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This work provides a unified theoretical treatment of the single and correlated double-electron emission from 
a general electronic system. Using Feshbach projection method, the states of interest are selected by the pro¬ 
jection operator; the Feshbach-Schur map determines the effective Hamiltonian and the optical potential for 
the emitted electrons. On the other hand, the nonequilibrium Green’s functions method is demonstrated to be 
a complementary approach and an explicit correspondence between both methods is established. For a self- 
contained exposition some results on single electron emission are re-derived using both formalisms. New in¬ 
sights and results are obtained for the correlated electron-pair emission: This includes the effective two-electron 
Hamiltonian, the explicit form of the Feshbach self-energy in terms of the many-body self-energies, and the di¬ 
agrammatic expansion of the two-particle current. As an illustration of the diagrammatic technique the process 
of the two-particle emission assisted by the excitation of plasmons is explicitly worked out. 
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I. INTRODUCTION 

Scattering experiments deliver the most detailed informa¬ 
tion on the structure of matter. For instance, the fully resolved 
spectra of an electron emitted from an electronic system upon 
photon or particle impact encode the spin and momentum- 
resolved spectral properties of the sample [1-5]. For direct in¬ 
formation on the two-particle properties the detection of a cor¬ 
related electron pair is necessary which is usually performed 
in a one-photon double-electron emission [4] or in a swift 
particle-impact double-electron emission experiment [6]. Cal¬ 
culations of the electron emission spectra from atomic and 
molecular systems [1,4, 7-9] as well as from condensed mat¬ 
ter [1-3] are done routinely. The underlying theories and tech¬ 
niques differ, however. The issue addressed here concerns the 
formulation of a unified and numerically accessible theoret¬ 
ical framework of single and double photoelectron emission 
(SPE and DPE) from finite and extended electronic systems. 
A method of choice for this purpose is the nonequilibrium 
Green’s functions (NEGE) approach [10-13]. In full general¬ 
ity the response function describing electron emission is more 
involved than the optical response which is related to time- 
ordered particle-hole (p*h) Green’s function (GE) for which 
well established approximations exist. Even for a single elec¬ 
tron emission the response function can only be defined on 
the Keldysh contour and after performing the calculations, the 
times are projected on the real observable times. The second 
complication is that for a fixed energy and momentum of the 
detected electron the sample maybe left in an excited state. A 
typical example is the plasmon satellites in core-level photoe¬ 
mission [14]. There, the target is left with one excited plas¬ 
mon [15]. The conservation of energy and momentum allows 
to focus on, e. g., the no-loss current. The response function 
is then determined by the product of two vertex functions and 
three single-particle Green’s functions [16]. If an approxima- 
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tion is made for one of the constituents, it has to be taken over 
consistently to the others. The notion of a conserving approx¬ 
imation is rooted in this requirement. 

Eirst theories of electron emission were empirical: E.g. for 
surfaces, following Berglund and Spicer [17] the photoemis¬ 
sion is regarded as a three stages process: excitation, transport 
to the surface (during this stage the particle may loose energy), 
and the transformation into a scattering (detector) state. In 
1970 Gerald D. Mahan wrote “we have not yet been able to de¬ 
rive a simple, time-ordered, correlation function which would 
serve as the starting point for a closed-loop type of calculation. 
That is, we have not yet found a ”Kubo formula for photoe¬ 
mission.”” [18]. Shortly thereafter Schaich and Ashcroft [19] 
and Langreth [20] employed a time-ordered formalism for the 
response function, and Caroli et al. [21] introduced the nowa¬ 
days standard NEGE formulation. The well-known Eermi 
Golden rule expression for the photocurrent 

Jp^2n f dedisp -s- w)(;if[r'|AA(e)AVp”^) 

%J —oo 

derives rigorously from the response-function formalism. In 
1985 Carl-Olof Almbladh obtained the following modifica¬ 
tions of the no-loss current; 

rti 
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xk\e-va),e)\Xp'’). 

In these formulas an interaction with an electromagnetic field 
of the frequency a> is assumed, ^ denotes the final scat¬ 
tering state with the momentum p and energy Sp, and A(e) 
is the spectral function. A(e H- u), e) is the so-called vertex 
function which, for noninteracting systems, reduces to the op¬ 
erator of the light-matter interaction A. In interacting systems 
it describes the screening of the optical field by the sample 
electrons and the accompanying polarization effects [22], 

The physics beyond no-loss has many facets. There are two 
prominent examples: the plasmon satellites [15, 23, 24] and 
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the Auger effect [25-28]. In both cases the system is left in an 
excited state that relaxes subsequently either due many-body 
effects or results in the emission of a secondary electron. It 
should be noted, however, that the borderline in such a clas¬ 
sification is blurred; one can consider the Auger effect as a 
two-step process, in which the decay is treated independently 
from the primary ionization or as the no-loss double photoe¬ 
mission [29]. The former point of view yields a description of 
the Auger effect in terms of an equilibrium two-hole Green’s 
function [26, 30, 31]. 

The goal here is to generalize the nonequilibrium approach 
as to treat single and double electron emission. We will mostly 
discuss processes related to the absorption of one photon. Par¬ 
ticle impact is discussed only in the optical limit as specihed 
in the Appendix A. In particular, this work provides a de¬ 
tailed discussion of DPE, a process that was experimentally 
realized for various systems [4, 32]. For a self-contained 
presentation we start by dehning observables and introduc¬ 
ing basic formulas solely based on the time-dependent per¬ 
turbation theory and the assumption of adiabatic switching of 
the light-matter interaction (Sec. II A). Already on this level 
one can reformulate these expressions in the Fermi golden 
rule form and demonstrate how the sudden approximation 
can be used to reduce the many-body to two-body description 
(Sec. IIB). Such reduction, however, neglects the energy loss 
of an emitted electron on its way to detector. These extrinsic 
losses are treated by means of the projection operator tech¬ 
nique (Sec. III). For single photoemission (SPF) this approach 
was established in works of Almbladh [16], Bardyszewski and 
Hedin [33], Fujikawa and Hedin [34], Hedin, Michiels and 
Inglesheld [35], and for DPF by Brand and Cederbaum [36]. 
The notion of the optical potential is central to this approach. 
While the case of elastic scattering was considered in a classi¬ 
cal work of Bell und Squires [37], the inelastic case, which 
is especially relevant for photoemission, is more involved 
and has a long history with a recent progress due to Ceder¬ 
baum [38, 39]. In Sec. IV we closely follow the derivation 
of Almbladh and extend the theory to the two-electron case. 
There are important differences as compared to the single¬ 
electron emission. Under some assumptions DPF is only pos¬ 
sible for interacting systems [40]. We demonstrate that the 
vertex function is the source of this electronic correlation ef¬ 
fect. Finally, we corroborate our hndings by performing a di¬ 
agrammatic expansion of the derived DPE response function 
in terms of Green’s function on the Keldysh contour (Sec. V). 
We consistently use atomic units. 


II. THE TWO-ELECTRON CURRENT 

For DPE from atomic and molecular systems [41, 42] a 
variety of very successful techniques, based on a full nu¬ 
merical solution or using approximate correlated scattering 
states of the few-body Schrodinger equation, were put for¬ 
ward. The wave-function-based methods and, consecutively, 
the scattering approach are less suitable for extended degen¬ 
erate fermionic systems. Such DPE experiments were hrst 
performed for Cu(OOl) and Ni(OOl) crystals [32] and mean¬ 


while for a variety of other samples. Here comes the response 
formalism into play: the expectation values of products of 
the creation and annihilation operators are computed over the 
ground state of a (many-body) system, and perturbative ex¬ 
pansions are evaluated with the help of Wick’s theorem. If 
the studied process can be regarded as a multi-step event, then 
the rate equations are often a very efficient tool. They can be 
derived either from the density matrix or from the NEGF for¬ 
malisms using some additional assumptions. For instance the 
generalized Kadanoff-Baym Ansatz has been used to derive 
the quantum master equations starting from NEGF approach 
to describe the transport in molecular systems [43]. 

Here we present a self-contained derivation of the two- 
particle current starting from the time-dependent perturbation 
theory. The resulting formula (Eq. (12)) is, however, less use¬ 
ful for practical applications because it requires (generally un¬ 
known) many-body states. One has either a choice to com¬ 
pletely neglect the target-ejected particles interaction which 
still might be relevant for higher energies (Sec. IIB), or, as 
will be demonstrated in the next section (III) to properly re¬ 
duce the formulations as to work with effective residual inter¬ 
actions (i.e. optical potentials). 

A. Basic definitions 

a. Hamiltonian: A system of interacting fermions is 
considered that has the Hamiltonian 

H - J"dxij/\x)h{x)>j/{x) 

where the field operator (0-') with argument x s (r, cr) anni¬ 
hilates (creates) a fermion in position r with spin cr. Needed 
below is the anti-symmetrized interaction 

y(xi,X2,X3,X4) = v(ri,r 2 )( 5 (X 2 - X3)6(xi - X4) 

- 6{x\ - X3,)6 (x2 - xa)). (2) 

One may wish also to change the basis for the representation 
of creation and annihilation operators via 

= '^{x\i)ci, (3) 

i 

where the sum runs over a complete set of one-particle states 
and we consistely skip ■' ■ on c, and c|. To study photoemis¬ 
sion we need to further classify the states according to their 
geometric character. A state will be called bound (0, e !B) if 
for any e > 0 there is a compact set B c such that for all 
times t the state remains in B\ 4>i\\ < where B‘^ is the 

complement of B, xb‘ denotes the corresponding characteris¬ 
tic function. Analogically for the scattering states (0k G C) we 
adopt the following dehnition; they are the vectors for which 

limj-^oo 2T J_j.\\XBe“^<f>k\\dt - 0 for all compact sets B c R^, 
i.e. they leave any bounded region. It is clear that !B ± C 
and according to the RAGE theorem [44] all the states from 
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the discrete (point) spectrum are bound, whereas the contin¬ 
uum states (absolutely continuous and singularly continuous) 
are the scattering states. Thus, parallels between the geomet¬ 
ric and the spectral classification allows us to use continuum 
and scattering, and point and bound terms interchangeably, 
although for the purpose of the present work the geometric 
classification is preferred. Finally we note that if our theory 
is to be applied to solids the use of localized Wannier func¬ 
tions [45] is preferred, at least for systems where their exis¬ 
tence can be proved [46]. 

We will use the letters {abed) for general orbitals, (/ jnm) 
for bound orbitals and bold-face letters for continuum states. 
In these notations: 

(4) 

ab abed 

^ ] ^ah^a^b ^ ^ abed^b^^ 

ab abed 

b. Initial state preparation: The above Hamiltonian de¬ 
termines the quantum state of the target (wave-function I'Pq) 
with corresponding energy Eq) in the remote past (f = -oo). 
When the system is perturbed by the interaction with external 
fields it evolves to a new state. As a typical mechanism we 
consider here the light-matter interaction 

y(f) = (Ae-“^‘ + K)e'‘^)e'^\ A = ^ AabcW (6) 

ah 

In this expression y(f) is adiabatically turned on allowing to 
introduce a typical interaction time ~ ( 2 ? 7 )“'. The form ( 6 ) 
permits generalizations: In Appendix A we consider the pro¬ 
cess of impact ionization caused a charged projectile particle 
(e.g. an electron) impinging on the target system. At high 
energy the projectile can be regarded as distinguishable from 
electrons of the system. This allows to average the projectile- 
target interaction over the projectile’s states and write the per¬ 
turbation in essentially the same form as in Eq. ( 6 ), i.e. as a 
single-particle operator. 

From the first-order time-dependent perturbation theory we 
obtain the approximate eigenstate of the full Hamilto¬ 

nian H + V{t) at time t -Q\ 

= I'Po) + lim--AIT'o). (7) 

Eq + 0 ) - H + irj 

Readers will immediately notice parallels of Eq. (7) with 
the scattering theory where the Mpller operators convert 
an eigenstate of H (the Hamiltionian of the target system) at 
t — +O 0 , into an eigenstate of H -H y(0) (the full Hamilto¬ 
nian) I'l'i’"^) = at time t = 0 (cf. Eqs. (14.66) of 

Joachain [47]). The scattering theory is required when elec¬ 
tromagnetic fields are quantized. For classical fields Eq. (7) 
follows from the first order expansion (in A) of the Mpller op¬ 
erator To emphasize the similarity we denote the state 
given by Eq. (7) as the scattering state. In what follows we 
omit the tilde which we used to denote its approximate char¬ 
acter. 


c. Observables: Assuming we know the quantum state 
of the target at f = 0 some observables can be computed. 
Since we are interested in photoemission these are the expec¬ 
tation values of the current operators. The safe way to in¬ 
troduce them is to use the continuity equation which is gauge- 
invariant. The one-electron current 7k is defined as the number 
of electrons Ak with a given momentum k outside the target 
divided by the effective interaction time ( 277 )“'. There is a 
detailed discussion [16] on why electrons in the sample give 
a negligible contribution to the current. Same arguments are 
valid for the two electron case. Thus, we analogically define 
the two-electron current as 


A,,k2 = lim277<Ak,Ak2 - (Jki.kjiVki)- 

77—>0 


( 8 ) 


In the expression above (and all subsequent derivations) we 
do not explicitly spell out the spin quantum numbers. The 
dependence on the spin can be recovered by substituting the 
continuum quantum numbers like k by kcr (likewise for bound 
indices). The second term excludes the one-electron current 
in the case when two momenta are equal. Eq. ( 8 ) gives access 
to the differential cross-section through the following relation: 


d^cr 

7ki(ik2 


a> 

7 


7ki,k2> 


(9) 


where I/uj is the photon flux density [48]. For the velocity 

gauge A = ^Ao ■ p, / = 377 , where Aq is the amplitude of 
the vector potential. Similar expressions can be given for the 
length gauge. 

The average in Eq. ( 8 ) is performed over the perturbed 
state (7): 


Jk„ki = lim 277<'1 'o|a1 -- — 

Eq + to — H — ITj 


.t 


ir 
I 


Eq + to — H + ir] 


-A|T'o>, (10) 


where we used the usual anti-commutation relations for the 
fermionic operators. The current is quadratic in A or linear in 
the number of absorbed photons. The first order in A gives the 
linear conductivity current and is of no interest here. [ 21 ] 

To derive the Fermi golden rule for DPE we insert a com¬ 
plete set of the (N - 2)-particle states and use the scattering 
theory to evaluate matrix elements of the type: 


Mi 


ki,k 2 ,;S 


= <'To|Ai 


1 




Eq + CO - H - ii] 


K)- 


We will generally use lower indices to distinguish quantum 
states and upper indices to indicate the charge of the system 
or the nature of the state (+), i.e., incoming or outgoing. For a 
scattering process with the following energy balance 

Ei = Eq + CO ^ Ef ^ Ski -H Skj + 


the Mpller operator 1 translates a wave-function in the re¬ 
mote future into a incoming (they are sometimes called in¬ 
verted FEED states [35]) scattering state at f = 0: 




-IT] 


Ef - H - iij 


.J J 

^kAk 2 


K)- 
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Following Almbladh [16] we obtain; 




resulting in the Fermi golden rule for DPE for an adiabatic 
switching of V(t): 


Ai,k2 = lim277^|Mki,k2,/?f 

P 

= 2nY^6{Ei - £/)|<'F«|A|'I'o>f. (12) 

P 


This is essentially an exact equation if strong field effects 
are neglected, i.e. if the first-order perturbation theory in field 
strength is adequate. Now we discuss some common approx¬ 
imations. In the sudden approximation the Mpller operator is 
set to the identity operator and it follows ss 
leading, e.g., to Eq. (1) of Napitu and Berakdar [49]. The sud¬ 
den approximation is broadly used to interpret the single pho¬ 
toemission. However, it is easy to construct an example when 
it completely fails: Consider photoemission from a system 
surrounded by a impenetrable potential barrier. Irrespective 
of the photon energy there will be zero current in the detector. 
Thus, it is extrinsic losses [35] that are missing in the sudden 
approximation. 


B. Sudden approximation 

In the sudden approximation for SPE it is possible to reduce 
the many-body description to a single-particle picture which 
also allows to approximately treat the Mpller operator and ac¬ 
commodate extrinsic losses. The central object in such an ap¬ 
proach are the Dyson orbitals [50]. The hole Dyson orbital is 
defined as an overlap of {N - 1) many-particle state with the 
N-particle initial state: 

0„(xi) ^ y/N J'd(x2 . ■ . Xiv)['P+(x2 ,..., x^v)]* 

X T'oCxi,... ,xm) = <T':|iA(xi)|'Po), (13) 

A rather extensive review of such overlap operators as well as 
the proof on the last ’’dressed in the fancy outfit of the occu¬ 
pation number formalism” identity can be found in Ref. [51]. 
Practical approaches for their computation are overviewed in 
Refs. [52, 53]. By introducing a similar two-hole Dyson or¬ 
bital. 

<l>f\xi,X2) = J"^(-^3 • • ■ Xw)[T'^^(X3, . . . , Xa,)]* 

X T'o(xi,... ,xm) = ^<'P| 1f(^i)0-(x2)|'Po>, (14) 

V2 

and neglecting the Mpller operator we obtain for the two- 
particle current (12): 

7k„k. - 2nJ]6(Ei - £y)|<kik2|A|0®)f, (15) 

P 


where |kik2) is asymptotic two-particle state, i.e. anti¬ 
symmetrized product of two plane-waves. The two-hole or¬ 
bital is anti-symmetric with respect to the interchange of par¬ 
ticle coordinates and in general has norm < 1. To derive (15) it 
is instructive to consider first a corresponding matrix element 
for SPE: 

Mk,, - F ^ . yA,,<>T:|ckclc,|>To). 

Now we have Ckcici|'Po) = bk,flCi|'Po) H- cic/,Ck|'Po) and it is 
time to make another very important assumption: 

Ckl'Po) ^ 0. (16) 

It is not valid in general, however, one can use the same ar¬ 
guments as Almbladh (see discussion around his Eq. (11)) to 
demonstrate that it gives a vanishing contribution. Eor homo¬ 
geneous electron gas this is even a generally valid statement. 
Besides allowing to compute the matrix elements the assump¬ 
tion (16) also justifies why terms resulting from the second- 
order perturbation theory give vanishing contributions to the 
current. 

In this way (see Appendix D) Mk^a = g_g^^.^ (k|A|0Q,) and 
7k = 2;r^b(£, -£/)|<k|A|0 „)f. 

a 

Eor DPE we analogically analyze the matrix element entering 
Eq. (11) and neglect terms with two holes at momenta kj and 
k2 (i.e. Ck^CkJ'Po) ~ 0) as compared to the terms with only 
one hole (Appendix D). Notice that for SPE we neglected one 
hole term as compared to zero hole contribution (cf. Eq. (16)). 

It is obvious that the sudden approximation is only valid for 
large momenta 2 and it is indifferent to the state in which the 
system is left in (the final double ionized state can be an ex¬ 
cited state). Thus, it is desirable to generate improved approx¬ 
imations to Eq. (12) by rewriting it in the two-particle form, 
but with an improved final state (such as Eq. (4) of Eominykh 
et al. [54] or Eq. (2) of Eominykh et al. [55]). 


III. EXTRINSIC EFFECTS 

A many-body target interacts with light such that certain 
number of electrons is emitted. Here, the fundamental ques¬ 
tion is whether it is legitimate to describe the process in such 
a way that only quantum numbers of ejected particles are con¬ 
sidered and remaining degrees of freedom are traced out, i.e. 
put into some effective interactions. The projection operator 
formalism is a general method to treat this kind of problems. 
In this section we introduce the basic concepts of this theory 
and demonstrate the reader that a deep connection with the 
nonequilibrium Green’s function formalism exist. We con¬ 
clude this rather mathematical section by considering two ex¬ 
amples. Based on these examples the Fermi golden rule is 
derived in the subsequent section. 
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A. Nonequilibrium Green’s fnnctions 

In the Keldysh formalism [13] the field operators evolve 
on the time-loop contour C shown in Fig. 1. Operators 
on the minus-hranch are ordered chronologically while op¬ 
erators on the plus-branch are ordered anti-chronologically. 
Letting z\ and zi be two contour-times, the Green’s func¬ 
tion G{x\Z\, X 2 Z 2 ) can be divided into different components 
G°‘^{x\ti,X 2 t 2 ) depending on the branch a,(3 - +f- to which 
Zi and Z 2 belong. As before, Xj, denote a composite coordinate 
comprising space and spin variables. For a = yS = - we have 
the time-ordered Green’s function 

G^^{x\ti,X2t2) = -i{T (17) 

In this expression the average {...) is taken over a given den¬ 
sity matrix p and T is the time-ordering operator. The sub¬ 
script “i/” attached to a general operator O signifies that that 
operator is in the Heisenberg picture 

dH{t)^€i{to,t)d€{{t,to), (18) 

where 'U(ti , 12 ) is the time-evolution operator and fo is an ar¬ 
bitrary initial time. Reversing the time arrow the G is con¬ 
verted into the anti-time-ordered Green’s function 

G**{xih,X2t2) = -i{T [^H{Xlh)ll/'^u{X2t2)\), (19) 

where T orders the operators anti-chronologically. Finally, 
choosing zi and Z 2 on different branches we have 

G^*(xituX 2 t 2 ) = i{!jr\{x 2 t 2 )^Hixiti)), ( 20 a) 

G*^{Xiti,X2t2) = -i{^H{x\tx)^'^jj{X2t2))- ( 20 b) 

The last two components are equivalently written as G * - 
G^ (lesser Green’s function) and G* = G^ (greater Green’s 
function), and describe the propagation of an added hole (G^) 
or particle (G^) in the medium. 

It is often convenient in addition to time ordered and anti¬ 
ordered functions to introduce the retarded and advanced com¬ 
ponents; 

G’^(xi, V 2 ; t) - 6(t) [G^(xi,x 2 ; t) - G^(xi,x 2 ; 0] > ( 21 a) 
G^(xi,x 2 ;f) = Q(-t) [G^(xi,x 2 ;f) - G^(xi,x 2 ;f)] ■( 21 b) 

In order to find their representation in frequency space we 
multiply the retarded GF by e ''' with ?; —> O-i- in order to 
enforce the convergence and compute the Fourier integral; 

G^(Xx,X2-,(jJ') = {^(Xl) -- —^\x2)) 

(jj + Eq - H + vq 

+ ^(X 2 )-— -Hx\)) ■ ( 22 ) 

Eq- CO - H -iTj 

t 

— 00 _ w 

- ) —► 


FIG. 1. The Keldysh time-loop contour C. The forward branch is 
denoted with a label while the backward branch is denoted by a 
“-H” label. 


Let further introduce (for general z 6 C) the particle-type and 
hole-type GF by 

G'p^(xi,X 2 ;z) = {^(xi)—^^\x 2 )), (23a) 

z- H 

G®(xi,X2;z) = {^\x2)—^il/(x\}). (23b) 

z- H 

From Eqs. (23) follows 

G’^^^(xi, X2; w) = G*-P^(xi, X2; Eq (jj± irf) 
-G®(xi,X2;£’o - o + \ t ]). 

Finally, let us present the equation of motion (EOM) for the 
retarded GE in the form; 


(u) \rj)G'^(x\,X 2 ', co) - d(xi - X 2 ) 

+ ([i^(xi),H] 


Eq + CO - H irj 


— ¥(X2)) 


<(A'TX 2 ) 


1 


Eq+ (jO- El + 12] 


-[ilf(x 2 ),H]). (24) 


The two-particle Green’s functions are much more diverse. 
However, we will only need those containing creation opera¬ 
tors with the same time argument and the same holds for anni¬ 
hilation operators. To specify the relative order of creation (or 
annihilation) operators infinitesimally small times are added. 
Because such Green’s functions depend on two times only, the 
same nomenclature as in the single-particle case can be used. 
Thus, we define 

G'PP^(xi,X2 ;Xi,X2 ;z) = {\j/(xx)^(X2)—^^\x2)^\xi)) , 

z-H 

G‘''''^(xi,x2;xi,x2;z) = {^\x2)^\xi) - -^(xi)ij/(x2)) ■ 

z-H 

They are the constituents of the retarded and advanced two- 
particle Green’s functions; 

/G^^^(xi,x2 ;xi,x2 ;<y) = G^pp^(xi , X 2 ; xi, X 2 ; Eq + 00 + it]) 
-G'‘'*'>(xi ,X 2 \x\,X 2 \Eq -oo + irf). 


For the retarded function the following equation of motion can 
be derived; 


(co ■¥ i77)G’^(xi, X 2 ; xi, X 2 ; w) 

- 6(Xi - Xi)G^(X 2,X2,0) - (5(Xi - X2)G^(X2 ,Xi,0) 

-H 6 (X2 - X2)G^(xi,Xi,0) - 5(X2 - Xi)G^(xi, X2, 0) 

- i{\^(x\)!l/(x2),H\ -—^l^■'■(x 2 )iA'^(xl)) 

^ Eq + CO - H + IT] 

- i{\j/Hx2)!lf\xi) -^ [t^(xi)iA(x2),i7]> . (25) 

Eq — CO — H — \rj 


B. Two projection operators 

In the previous section we have seen that relevant types of 
Green’s functions can be written in the form of a resolvent 
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((z - z e C. To be more specific about the state over 

which the averaging is performed we select from all possi¬ 
ble states of the target and emitted particles the relevant ones 
for the effect of interest by employing projection operators. 
In the following we consistenly skip ■' ■ when writing these 
operators and use 1 to denote the identity operator. Hence 
P + Q - 1 are two complementary projection operators with 
the idempotence {P^ = P, - Q) as their defining property 
and the basis formula for computing resolvents 


1 


z- PI z- Pip- £/>(z) 


1 H- PHQ- 


1 


Z-Hn 


(26) 


where Hp - PHP, Hq - QHQ, and the self-energy operator 
is defined as: 


tp(E) = PHQ - ^-^QHP 

E - Hq 


(27) 


The map Fp \ H ^ Hp + f,p(E) is called the Feshbach-Schur 
map, it relates the eigenvalue problem on the full Hilbert space 
and to that on its subspace. We summarize relevant matrix 
identities in Appendix C. Due to the presence of the bath 
Hamiltonian Hq in Eq. (27) this definition cannot be used 
for practical computation of the self-energy. Fortunately, a 
connection with the many-body perturbation theory (MBPT) 
exists [56, 57]. If, for example, starting from the A-particIe 
Schrodinger equation //|'I'o) = Eol'To) we use a projector 


^a\^) 


where n(r) is the hole-density of ionized state a, i.e. 
h(r) s {'l'+|0'(r)i^'(r)|'l'^), the eigenvalue problem on the P- 
subspace (C3) «'I'+|i/^(r)(i/p -H tp{E) - £/f)P|'I'o) = 0) is the 
Lipmann-Schwinger equation for the hole Dyson orbital (13). 
Notice that Hp contains the electrostatic and exchange part of 
self-energy, whereas t.p{E) —> 0 for E —> +oo. Similarly, in 
1959 Bell and Squires [37] considered a one-body potential 
for the scattering of a particle incident on a complex (many- 
body) target. They demonstrated that this optical potential is 
exactly given by the sum of all proper linked diagrams, i.e. 
many-body self-energy in the time-ordered formulation. In 
fact, their Eq. (7) directly corresponds to Eq. (C3) when P is 
a projection yielding a particle Dyson orbital. 

In order to study single and double photoemission we intro¬ 
duce two special projection operators. The main goal of this 
section is to establish an equivalence between the abstractly 
defined self-energy (Eq. (27)) and the self-energy of the many- 
body perturbation theory. We consider the expression appear¬ 
ing in the first line of Eq. (26) i. e. resolvents of the type 

1 1 
P - -P = P -^^- P- 

z-H z- Hp- 2p(z) 


We will demonstrate that the formalism of nonequilibrium 
Green’s functions is easily paralleled with the Feshbach pro¬ 
jection algebra (FPA). The basic relation for the subsequent 


derivations are the operator identities 

(A - B)-^ = -H A-^B(A - B)-^ , (28a) 

(A - B)-^ = A^' +(A- By^BA-^ . (28b) 

We will show below that with 

A = z- PHP = z-Hp, (29a) 

B = QHP + PHQ + QHQ , (29b) 

the operator identity (28) has a structure of the Dyson equation 
for certain Green’s functions. 

For SPE we consider the projection operator 

p.(30) 

k 


where the sum runs over scattering states. It is common to 
select these single-particle states |i,cJk) to be eigenfunctions of 
some reference Hamiltonian with proper boundary conditions. 
We request that I'P^) is a completely bound remainder of the 
ionization event and does not emit a second electron at a later 
stage (Auger electrons are a typical example for these kind of 
processes). There are many equivalent ways to impose this 
restriction, for instance we will assume 


Ckl'T:) = 0, (31) 

i.e., implying |'P^) is a vacuum state for photoelectrons. From 
the assumption follows the indempotency (P^ = Pa, see Ap¬ 
pendix D for proof) and, thus. Pa represents a true projection 
operator. The application of Pa restricts the possible processes 
which might occur upon excitation to the definite emission 
of one photoelectron, whereas the ionized system is left in 
a (possibly excited) bound state I'T^). From the assumption 
Eq. (31) follows another restriction: 

lim - lim y{x\i}ci(t)\%} 

r—ioc r—ioo t * 

i 

+ lim y<x|k>Ck(f)l'P:) = 0, (32) 

r—ioc / A 

k 

where the first term is equal to zero because each bound state 
(/) is necessarily given by a square integrable function (con¬ 
verse is not true). In the following we will use another conse¬ 
quence of assumptions Eqs. (31, 32): 

GL(") = 0, = 0, (33) 

lim G'^(xiti,X 2 t 2 ) - lim G'^(x 2 t 2 ,xiti) - 0. (34) 

rj—»oo ri—»oo 

The projection operator for DPE we define as 

(35) 

pp' 

Here, is the doubly-ionized reference state, to which two 
photoelectrons with continuum quantum numbers p and p' are 
added. We can easily show the indempotency of the projection 
operator (35) if we require, similar to Eq. (31), 

Cfl'i'p = 0 . (36) 
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C. Example of SPE 

d. Equation of motion (EOM): As a starting point let us 
use the following operator identity which can be derived from 
Eq. (28a) or verified by direct computation 

(Z - E^JPa-^Pa ^Pa+ Pa(H - E^J-^P^ . 
z-H z-H 

With the definition of the SPE projection operator P^ in 
Eq. (30), we find 

Z pq Z n 

pq 

where we applied the definition of the particle-type GE 
Eq. (23a). Note that the GE is defined for a particular sub¬ 
space spanned by the operator Pa and should therefore always 
be understood as the GE associated to |'P^). Eor brevity, how¬ 
ever, we omit labelling GE by a. 

Using these notations the operator identity reads 

(z - e:) 4i'i':)GX\z)<'i':icq = ^ 4i'p:)<'p:ick 

pq k 

pq Z M 

With the help of our assumption Eq. (31) we can now remove 
the sum by applying ('P^kp' from the left and from 

the right as Eq. (31) implies ('Pakp/CpI'P^) = 5pp'. Eurther- 
more, we note that f¥^\cp(H - Ef) - ('I'^|[cp,i/] because of 
H\'i‘a) - Hence, we obtain 

(z - £:)GX\z) = <5pq + [cp, h] . (37) 

As stated above, we can think of |'P+) as a vacuum state for 
free particles (cf. Eq. (31)). The hole-type GE is identically 
zero. Therefore, 

G;i;\£: + w + i77) = G^q(cu), 

Substituting z = Ea + a> + iq in Eq. (37) we realize its equiva¬ 
lence to Eq. (24). In other words, by applying the EPA we can 
derive EOM for the retarded Green’s function. 

e. Effective Hamiltonian: In Eq. (28a) plays the 
role of the reference Green’s function. Correspondingly, 
PHP is the effective Hamiltonian. Using the standard anti¬ 
commutation algebra and the assumption (31), we find 

('P^lcpi/cJlT':) = £>pq + <'P:| [cp, h] cl\%) 

~ Pa^ptl ^pq ^ (pp«mq ~ P«pmq) 

nm 

~ Pa^pti ^pq’ (38) 


i.e. it consists of the total energy of the ionized system and the 
Hartree-Eock Hamiltonian for continuum states. The latter is 
computed with the density matrix of the target: 

fpq “ ^pq "t" ^ ] ’T’p/jqm(^n^m)- (39) 

(nm)£&- 

Let h be an operator acting on the subspace of continuum 
states with matrix elements given by Eq. (38). Its resolvent 

^p’’q(2) = <'I':kp^4|T':) (40) 

z - h 

relates to the reference retarded GE as gpq(w) = gpq^E^ +a) + 
iq). 

f. Self-energy and the Dyson equation: The second cor¬ 
relator in the EOM (37) amounts to 

Z n ^ Z n 

+ ^ ^ l^pnabCPal^l^aCb ^‘^ql'^^o;) • 

n ab Z — H 

With Eq. (29) inserted into the identity Eq. (28a) we apply Pa 
from left and right, use the same trick to multiply with suitable 
states from left and right, and find 

Gpq(z) = glfq^z) - (Z) 

kk' 

k a 

+ Xi Xi X . 

k « afc Z H 

With z = Ea+(x> + iq Eq. (41) has a structure of a Dyson 
equation for the retarded Green’s function in the subspace of 
continuum states: 

G^q(^) = 4q(") + Z 4k(‘^)^L(")G«^(w). (42) 

kG 

The second sum runs over the full set of orbitals (bound and 
continuum). This is the most general form and without addi¬ 
tional analysis it cannot be reduced to the Dyson equation with 
the self-energy from the projection formalism (cf. Eq. (27)). 
Let us compare Eq. (41) and Eq. (42). At first we notice that 
Eq. (39) defines the reference Hamiltonian only on the sub¬ 
space of scattering states. We might extend the definition and 
request, for instance, that all the basis functions (bound and 
scattering) are the eigenstates of the reference Hamiltonian. 
This implies fpq = ep5pq and f„q = 0. Thus, mean-field terms 
of the Hartree-Eock Hamiltonian are then cancelled by the fre¬ 
quency independent part of the last correlator in Eq. (41). In 
the case when the reference Hamiltonian is not diagonal in 
the chosen basis the embedding self-energy terms addition¬ 
ally appear. In the simplest case (no interaction), they can 
be written as Epq(z) = f^mn tpng^nlizfnq- Let US now assume 
that the single-particle basis is such that no embedding self¬ 
energy appear. What would be the diagrammatic structure of 
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for the retarded GF from the many-body perturbation theory. 
The difference between them is the domain where the self¬ 
energies are defined; the Feshbach self-energy operates on 
the continuum sector only, whereas many-body perturbation 
theory does not impose such a restriction. By writing another 
Dyson equation (43) in the bound-continuum sector we can fi¬ 
nally obtain the Dyson equation with an effective self-energy 
in the continuum-continuum sector. This self-energy is an 
exact counterpart of the Feshbach self-energy (27). To the 
best of our knowledge it is the first explicit example of such 
correspondence. Critical for our derivation was the choice of 
the single-particle basis. We have demonstrated that it is the 
projection operator that determines the effective Hamiltonian, 
and if the basis is such that the Hamiltonian is diagonal the 
embedding self-energy vanishes and one arrives at Eq. (44). 
No further assumptions have been made and Eq. (44) is so far 
exact. 


FIG. 2. (a) Example of self-energy diagram that mixes bound and 
continuum states and is the building block of the second term in 
brackets in Eq. (44); (b) Mean-field Hartree contribution to the ef¬ 
fective Hamiltonian Eq. (38); (c) A typical contribution to the elec¬ 
tron self-energy in continuum-continuum sector in the case when the 
photoelectron is completely screened in the sample. 


the self-energy (27)? Erom the Dyson equation in the bound- 
continuum sector 


mk mn 

(43) 

we determine the Green’s function in this sector (Gbc) and 
substitute in Eq. (42): 


Gpq(w) = 4q'^w)-H^g^^(m) 

kk' 


"h ^ch 


8b 


1 - gb^bb 


■^bc 


kk' 


^k'q’ 


(44) 


where for brevity the subscripts b and c denote the bound and 
the continuum sectors. Expression in square brakets (Eq. (44)) 
can now be compared with the self-energy from the projection 
formalism (27). Notice, that the reference Green’s function 
was assumed to be diagonal, i.e. gb = ghb and gbc = 0. 

g. Dominant scattering mechanisms: Let us recapitulate 
what led us to Eq. (44). We have chosen a projection opera¬ 
tor in the form (30). This specifies the state of a system af¬ 
ter the photoionization as containing one photoelectron in the 
scattering state plus the bound ionized target. Next, we ob¬ 
tained an effective Hamiltonian (38) acting on the P subspace 
and used it to define the reference Green’s function (40). We 
want to understand what is the diagrammatic content of the 
Eeshbach self-energy (27). It is not possible to use this equa¬ 
tion directly because it involves the effective Hamiltonian on 
the complementary Q-subspace. However, it is possible to 
use another matrix identity (28a) and to formulate the Dyson 
equation for the full Green’s function in the P subspace (41) 
avoiding the use of the QHQ resolvent. This equation can 
be put in a direct correspondence with the Dyson equation 


Let us analyze the meaning of different terms of the 
photoelectron self-energy (Eig. 2). As discussed in details 
by Bardyszewski and Hedin [33], Almbladh [16] and Eu- 
jikawa and Hedin [34] scattering states vanish in the sam¬ 
ple (damped) represent the real photoelectron states more pre¬ 
cisely. One can derive explicitly the residual interaction that 
they experience. The reasoning is easier to perform in real 
space where the Coulomb interaction depends on two coor¬ 
dinates only (cf. Eq. (2)) as opposite to the Coulomb matrix 
elements which are four index quantities. Since the scattering 
states are damped in the sample, there are only two nonvan¬ 
ishing Green’s functions G™ and Gyv operating exclusively 
in the inner (v), outer (V) spaces, respectively. The Green’s 
function starting in the sample and ending outside of it (Gy,.) 
and the reverse (Gi,v) vanish. We can rewrite Eq. (44) in these 
new notations, however, it is not even necessary as it amounts 
to mere replacement 1? —> v and c ^ V. What has changed 
is the interaction lines in the diagrammatic expansion of the 
self-energy. They can connect v and V domains and generate 
therefore nonzero contributions. It is easy to see, however, 
that the second self-energy term vanishes: a diagrammatic ex¬ 
pansion of Si.y necessarily contains at least one g,,y line which 
is zero according to our assumption. Thus, only Eyy needs to 
be analyzed. By explicitly forbidding the particle exchange 
with the sample we arrived exactly at the case of elastic elec¬ 
tron scattering considered in the seminal paper of Bell and 
Squires [37]. We will see below that the structure of Syy is 
quite general and appears in the diagrammatic consideration 
of other processes, remarkably, in the parquet diagram treat¬ 
ment of the Eermi edge singularities [58]. There, however, 
a similar diagrammatic expansion arises due to the specific 
choice of the interaction between the deep hole (labeled by 
m) and the conduction electrons: Hi - Vkk'C^Ck'GjjCm. 
In contrast to their work, what induces a special structure of 
diagrams for Syy is not a specific form of the interaction ma¬ 
trix elements, but rather the absence of the off-diagonal blocks 
in g. It is easy to construct the electron self-energy fulfilling 
these restrictions: it consists of one open photoelectron line 
(depicted as solid line on Eig. 2 and a number of closed bound 
electron loops (depicted as dashed lines). Because of the re¬ 
striction (33) there are no photoelectron loops. 
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The topic of the present section is quite extensive and such 
an aspect as the Lehmann representation of the Green’s func¬ 
tions mentioned here was completely left out of our discus¬ 
sion. This is, however, very relevant for the treatment of finite 
systems, with important recent progress, e.g., [59]. 


final results read as 

(Z - = 5 pq( 5 p'q' - ^pq'^p'q 

+ [cpCp,, . (45) 


D. Example of DPE 


h. Equation of motion: The derivation for the two- 
particle case goes along the same lines. We insert the defi¬ 
nition of the projection operator (Eq. 35) in the identity 

(z - Ef)Pp^Pp = + Pp{H - , 


replace {'i‘f\cp,Cp{H-Ef) = <'T2+|[cp-Cp, H], and as for SPE 
compute the matrix elements of the whole expression. The 


The prefactor 1 /4 originating from the product of two projec¬ 
tion operators is cancelled because of the symmetries of the 
particle-particle GE and of the second term on the right-hand 
side of Eq. (45): 


/--(pp) 

p'pq'q 


(z)^G 


(PP) 

pp'qq' 


(z) = -G 


(PP) 

p'pqq' 


(z) - -G 


(pp) 

pp'q'q 


(z) ■ 


(46) 


Inserting z = + lo + \q shows the equivalence of Eq. (45) 

to the equation of motion (25). 

i. Effective two-particle Elamiltonian: Analogically to 
the SPE case we consider the Feshbach-projected Hamilto¬ 
nian in the subspace defined by Pp and describing two elec¬ 
trons including their interaction and their mean-field interac¬ 
tion with the ionized system: 


QPP |Cp'Cpi/CqCq,|'P^ ) — Ep (5pq(5p'q' <5pq'^p'q) + I [‘^P'‘^P’‘^q'^q'1*1/? )> 

where the last term can be expressed as follows 


(47) 


h'Cp,//] clcl,\%^} = fpqVq- 


-t- tp'q'^pq ^pq'^p'q tp'qdpq' 


+ Pi 


pp'qq' ^pp'q'q 


+ ['^Vnab{%^\clcp'CaCbc\c^^\'i‘4 - Pp'™i(T'|+|c^CpCaCicJc^, |'P|+)] . 

n ab 


(48) 


The first correlator in the square brackets evaluates in terms 
of the density matrix with respect to |'P^^) with bound state 
indices to: 


l^pnqm^p'q' ^pnq'm^p'qjitt^ttm} ■ 


Here we have written it in terms of the matrix elements of 
the anti-symmetrized Coulomb interaction (2) 'Vabcd = Vahcd - 
Vahdc- Similarly, the second correlator is obtained from this 
expression by the index exchange p <-» p'. The effective two- 
particle Hamiltonian (47) is so expressible as a Hartree-Fock 
Hamiltonian (39) for two independent electrons plus the inter¬ 
action (Fig. 3): 


h 


p'pq'q 


py (^pq<^p'q' ^pq't^p'q) + (fpq^p'q' + 7p'q'<^pq) 
(7pq'(5p'q + fp'q(Jpq'^ "I- 'T’p'pq'q- (49) 


P P' 



^1 Is 

q q' 


FIG. 3. Interaction between the photoelectrons incorporated in the 
effective Hamiltonian Eq. (49). Dashed lines denote bare bound 
state propagators. Dots denote the anti-symmetrized Coulomb in¬ 
teraction (2). 


and define the reference two-particle GF 


j. Kernel and Dyson equation: We return to the ma¬ 
trix identity (28a) and insert the splitting (29) with P - Pp 
(Eq. (35)): 


z-H 




z - h 


Pp + Pp 


1 


zPpH- 


1 


-Pa 


z-h z - El 
1 . 1 


z-h 


:PpHPp 


Z-H 


zPp^ 


^Eqq'^^^ = ('EjiCpCp, 


Z-h 




(50) 


Invoking again the symmetries (46), which also hold for the 
reference GF, and applying the same states from left and right. 
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we obtain 

(z) = (z) + V (z) 

pp'qq'^^'^ ‘^pp'qq'^^^ *pp'kk'''^^ 

kk' 

- I ■ (51) 

nn' 

It is instructive to divide the kernel entering the equation of 
motion (second line of Eq. (45)) or the Dyson equation (sec¬ 
ond line of Eq. (51)) into the terms containing higher corre¬ 
lation functions and those expressible in terms of two-particle 
GEs: 

[ckCksi/] - rkk-qq'fe) 

+ Z - 'k^^XqXz)) + X ^kk'..GX,(z). 

b ab 

The latter gives rise to the particle-particle embedding self¬ 
energy. We can now formally introduce the correlated 
frequency-dependent and the static kernels: 


given as a product of fully-interacting single-particle GEs, 
but rather is the full two-particle GE — the resolvent of the 
effective Hamiltonian (47) which includes the full electron- 
electron repulsion and the mean-held contribution from the 
ionized system. 

IV. FERMI GOLDEN RULE 
A. Single photoemission 

SPE was treated by several authors. We recapitulate the 
main points. The total observed current is proportional to the 
expectation value of the electron number operator iVk = c]^Ck. 
Out of all possible hnal states of the target we discard all un¬ 
bound states, i.e. Ckl'P^) = 0 and choose only those relevant 
for a specihc experiment. Let be a corresponding distribu¬ 
tion function. Eor instance when the target is left in the ground 
state we can set To = 1 and Aa = 0 for all excited states. Mod- 
ihed particle number operator for this process reads: 

Nk = ^T„4|T'+)<'I'+|ck = ^T,P,4ckP,. 


rkk-qq-(z) - X K-nnXz) + ^^k“k'„„']Gl’:!,qXz) 

nn' 

n ah Z H 

- Z Z '>^’nah{'Vl^\cWcaCb-^C^^C^^'^f) . (52) 

n ab Z n 

The static part is exactly cancelled by the density-dependent 
part of the effective Hamiltonian: 


*^kk'qq' “ ^ \ {^n^m}\^knqm^k'q' ^k'wq'm^kq 
nm 

~ k«q'm <^k'q ^k'wqm<^kq'j- (53) 


The embedding self-energy originates from the kernel as well 
as from the effective Hamiltonian (49): 


x~*(pp) / \ 

/ j ^kk'nn'^nn'qq ''^^' 
nn' 

= Z 

m 

+ ^ ykk'afcG|j^qq,(z) - ^ ykk'pp'G| 


(PP) 17') fzf'l 

kmqqA^.’ 'km'-fk'mqqA^ty 


(PP) 

pp'qq' 


(Z). (54) 


ab 


PP' 


The same expression can be obtained from the Langreth ap¬ 
proach starting from the Wigner distribution function [20]. 
Let now the SPE current be the expectation value of this op¬ 
erator 


7k = lim27,yT,<T'o|A' 


1 


7]-^0 


Eq + CO - H - ir/ 


P a 


X--A|T'o). 

Eq + o) — H + irj 


(56) 


We only consider the case 


1 


Ei-H + iT] Ei -Hp- t^;\Ei) 


(57) 


where we neglect the off-diagonal term in Eq. (26) and de- 
hne £^V<^) = £p(w + irj). We omit the subscript a where it 
does not cause a confusion. A simple calculation leads to the 
modihed matrix element 


Mk. = ('P+Ick- 


1 


Ei -Hp- t-;\Ei) 


-P.A|T'o>. 


(58) 


Using the same assumption for the computation of the matrix 
element of A, ('I'^|cpA|'Po) = {p|A|0Q,) and the dehnition of 
the Green’s function on the Pa subspace: 


With the results (49, 52, 53, 54) we can cast the Dyson 
Eq. (51) in the hnal form 

^pp'qq'^*^) ~ .?pp'qq'(*^) "*■ .?pp'kk'(^) 

kk' nn' 

X GLqq'(«) ■ (55) 

Eq. (55) has a form of the Dyson equation for the two-particle 
Green’s function, however, the reference GF gpp/qqX^) is not 


--(p) 

^pk,a 


{(o + ea± irf) = ('P^kp 


1 


Ei - Hp ■ 


■t^p\Ei) 


cliy) 


we obtain for the current 


7k = Yimlrj'^Aa ^<0„|A^|p)Gp'’^_^(w H- - irj) 
a pq 

X G^kq,ff(" + Sa + in){^\M‘Pa), (59) 
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where Eq, = Eq-E^. As shown in Appendix B we can express 
the particle Green’s functions in terms of Mpller operators 

+ s„ - iif) = _ -(60a) 


pk,Q''' 


OJ + Sa- S\^ - IT] 


■^(P) 


{u + Sa + irf) - 


1 


"“l’® ■' (0 + Sa 

This finally leads to the current 


Sk + 


-CfMlq)- (60b) 


4 = 27r V 4Crl 


(-) 


,)6(ci) + - ek)<<^a|A%ki)' 


A standard definition of the spectral function entails to 
A{O^Yj\4,,)6{^-ea){cl>a\. 

a 

Therefore, we can recast the expression for the current in a 
more familiar response form 


Using assumption (57) Eq. (62) can be written in the Fermi 
golden rule form with a modified matrix element 


^kik2,;S ~ ('^j8 ^k2t^ki 


1 


Ei -Hp- tf{Ei) 




Using the matrix elements of A, ('I'|^|cqCpA|'l'o) = {pq|A|0®) 
(cf. Eq. (D8)), and the properties of the two-particle Green’s 
functions (Appendix B) 






- 4 , 4 . 




( 2 ) 

w -H - Ek, - ek2 ± IV 


<pqi<L,P^ (63) 


4 = 27r r d(6(aj + ^ - Ek)(;rt.il^(OA'^lft4- 

%J —oo 

where the tilde denotes a spectral function with restrictions 
imposed by the weighting factors Aa and fi is the chemical 
potential, or in the Fermi golden rule form: 

7k = 2n AgSio) -H Eq. - ek)|(Tk4l^l^“)r' 

a 

The major distinction from other approaches is that both, ini¬ 
tial and final states are dependent on the final state of the target 
a. Formally, '■6® incoming scattering state of an elec¬ 

tron in the optical potential of the ionized target in the state 
I'T^). Notice that the current has been obtained using the ap¬ 
proximation (57). Exact calculation leads to the appearance 
of the vertex functions that describe a screening of the optical 
field by the electrons of the target [16]. We will stop on this 
point when treating DPE process. 

B. Double photoemission 

The total observed current is given in terms of the expecta¬ 
tion value of the electron number operators iVkjkj = A^k:A^k2 - 
4|,k2^k,, viz. Eq. (8). Out of all possible final states of the 
target we discard all unbound states, i.e. Ckl'T^^) = 0 and in¬ 
troduce weights 4 selecting the relevant ones. The modified 
observable reads: 

Nk2k2 = X'^/^4.42l'*'^"><'*'4l^k2Ck, 

= Xi ^^4i42Ck2Cki4. (61) 

p 

This allows us to improve upon Eq. (15): 


we finally obtain for Eq. (10) 

Aukz =27 t d^5(u + ^- Eki - ek2) 

\J —OO 

x<<L/^l^®(^)A'l<k2,4’ (64) 

where - max^(£o - EJ^) is the negative of second ioniza¬ 
tion potential, |i/'k 42 / 3 ^ *^62 incoming damped two-electron 

scattering state in the optical potential of doubly ionized target 
and A^^i(^) is the two-particle spectral function, which can be 
written in terms of two-hole Dyson orbitals: 

(65) 

P 

with E® - Eq- 

Notice that the current has been obtained using the approx¬ 
imation (57). Exact calculation leads to the appearance of 
the vertex functions resulting from 2ysA|'l'o) and describing a 
screening of the optical field by the electrons of the target [16]. 

In the valence shell the DPE mechanism is typically due 
to ground state electron correlation, i.e. due to the correlated 
two-particle spectral function entering (64). In contrast, when 
core electrons are involved a dominant mechanism for DPE 
is due to the final state relaxation (so called shake-off). Mul¬ 
tiple stages are then described by introducing corresponding 
projection operators for each intermediate stage. In the fol¬ 
lowing, we focus on the diagrammatic approach because it 
allows us to treat all these effects on equal footing. 


V. DIAGRAMMATIC APPROACH 


7ki,k2 = lim2;7^4<T'o|A 




1 


Ppclci^^k^C^,Pp 


Eo + to - H - irj 
1 


Eq + to — H + irj 


-A|'Po>, (62) 


Treatment of the off-diagonal part of the Hamiltonian resol¬ 
vent is the main difficulty of the Feshbach projection algebra. 
It is even more aggravated in the two-particle case. The dia¬ 
grammatic technique provides a natural and practical solution 
to this problem. 
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FIG. 4. Second order diagrams (in bare Coulomb interaction) repre¬ 
senting the DPE process. The dots labeled ki and k 2 correspond to 
the scattering state of two electrons observed in a coincidence mea¬ 
surement by the detector. Notice that not all combinations of pluses 
and minuses are possible because Coulomb interaction can only con¬ 
nect vertices on the same branch of the Keldysh contour, (a) Diagram 
vanishes according to the assumption (33) for dressed GFs. (b) Dia¬ 
gram vanishes because it contains an isolated island of minuses, (c) 
and (d) are the lowest order nonzero diagrams. The remaining two 
are obtained by permuting ki and k 2 . 

A. Derivation 

Eq. (10) when transformed to the time domain gives rise to 
the following ground state correlator: 

= <'I'ol4(f)c«(f)4^(0)4^(0)ck,(0)ck,(0) 

X 4(f')Q(t')l'I'o), (66) 

where the field operators are in the Heisenberg representation 
and f, f 6 (-cx),0] are physical times. For clarity, we omitted 
the indices in the notation of the correlator. It can be eval¬ 
uated diagrammatically by adiabatically switching on the in¬ 
teraction in the remote past, i. e. Hg - Hq + Now 

the average is performed over the noninteracting ground state 
|<I)o) and the times ^ lis on forward, backward branches 
of Keldysh contour y (Fig. 1), respectively: 

Z(f,f') = 

X <(0)42(0)ck2(0)ck,(0)4(fl)crf(f:))|<I)o>. (67) 

T here is the usual contour ordering operator [13] with the 
order relation x. i/a is such that it is equal to the Hamiltonian 
of noninteracting system Ho in the remote past and is iden¬ 
tical to // at f = 0. Notice that it is different from adiabatic 
switching on of the electromagnetic field in Eq. (7). |<l)o) is the 


ground state of Hq. Using Wick’s theorem we can contract the 
product of field operators in order to express the correlator in 
terms of products of single-particle Green’s functions. Zeroth 
order obviously yields four fermionic lines. However, if we 
use the same assumption as in Sec. IV B any zeroth order dia¬ 
gram vanishes. This is easy to understand by comparing with 
SPE case. There, no-zero contributions are coming from the 
following contraction: 


This is the only combination that results in greater GFs when 
one of the arguments is a scattering state (and is compatible 
with (33)). In particular, the above contraction equals to 

In DPE two creation operators with continuum state indices 
need to be contracted with two annihilation operators on the 
positive track. However, there is only one such operator. 
Hence, 0th order in interaction is zero. The argument that ex¬ 
cludes the first order diagram is slightly different and is based 
on the fact that bare interaction is instantaneous, i.e. corre¬ 
sponding time-arguments necessarily lie on the same, positive 
or negative, track. 

Second order nonvanishing contributions contain products 
of two Coulomb interaction operators (e.g. at contour times 
t+ and f_) and already a familiar product of six operators as 
in Eq. (67). From all possible contractions (they yield eight 
fermionic lines) we have to exclude many terms. Some of 
them immediately vanish because of the assumption (33) for 
noninteracting GF. Others, represent the Hartree-Fock renor¬ 
malization of two fermionic lines and likewise vanish because 
of the same assumption for the full fermionic propagators. 
Fig. 4 (a). Then, there are diagrams (Fig. 4 (b)) containing iso¬ 
lated islands of pluses and minuses which also vanish because 
otherwise the two-particle current cannot be written in the 
Fermi Golden rule form [60, 61]. Finally, there are only four 
(times two for exchange) nonzero diagrams. Two of them are 
depicted at Fig. 4 (c,d). 

It is clear now how more general diagrams for the two- 
electron current can be constructed: i) One replaces all bare 
fermionic propagators and interaction lines with the dressed 
ones; ii) Each pair of parallel fermionic lines are replaced by 
the corresponding two-particle propagator. Fig. 5 (a). In do¬ 
ing so one obtains, in principle, diagrams given by Fig. (lb) 
of Fominykh et al. [54] with a small correction that zeroth 
and the first-order two-particle GF should be excluded from 
the vertical track; iii) Next class of the diagrams are those that 
describe the screening of the optical field. Fig. 5 (b); iv) Pro¬ 
cesses involving intrinsic or extrinsic losses are given by the 
diagrams with interaction lines connecting points on different 
tracks, i. e. “-H “-H 0”, 0”. They cannot be obtained 

by the renormalization of fermionic or bosonic propagators, 
one such example shown at Fig. 5 (c) reveals a process with 
extrinsic losses. 

Finally, we give a description of a general diagram for a 
photoemission process. Examining SPE and DPE diagrams 
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FIG. 5. (a) Diagram for the two-particle current involving dressed 
two-particle propagators, (b) Simplest diagram where the optical 
field is screened, (c) Example of a diagram describing external 
losses. Thick wavy line denotes the screened Coulomb interaction, 
(d) Generic diagram for the two-particle current. 


we see that all of them are constructed from the common an¬ 
cestor: the density-density response function s + having 
a form of two islands with time arguments belonging to either 
forward or backward tracks of the Keldysh contour. Now we 
introduce detectors (shown as black squares at Fig. 5 (d) mea¬ 
suring 7ki,k2- As explained before i) the lesser GF with one 
of the indices being a continuum state vanishes because of 
the assumptions (31,32); and ii) observation is made at the 
rightmost point of the contour (i. e. at f_ = f+ = 0 in our no¬ 
tations), thus, each detector measuring particle numbers Nk, 
is connected to two greater GF. In view of this, the detectors 
“lie” on the fermionic lines flowing from the (forward 
track) to “H-” (backward track) islands. Each response func¬ 
tion constructed in this way has an important property that it 
can be represented in the Fermi Golden rule form, such con¬ 
struction obviously generalizes to an arbitrary number (n) of 
emitted particles. Simple counting shows that these processes 
are of at least 2(« - 1) order in the Coulomb interaction. 

The diagram in Fig. 5 (d) is a generic one describing all the 
DPE processes including the ones with losses such as shown 
at Eig. 5 (c). One can go a step further and give a prescription 
for classes of lossless diagrams. A detailed analysis of this 
particular situation is possible and will be done elsewhere. 
Here, we mention without a derivation that such diagrams 
can be split into the scattering part (the two-particle propa¬ 
gators can be written in terms of the scattering states llAk^k^^), 
cf. Eq. (63)) and the spectral part (containing the two-particle 
spectral function, Eq. (65)). 



FIG. 6. Diagrams for the plasmon assisted photoemission. SPE 
setup: only the primary (a), secondary electron (b) is observed, the 
fate of another electron is not specified, (c) DPE setup: both, primary 
and secondary electrons are observed in coincidence. 


B. Example of plasmon assisted DPE 

As an example we consider the processes depicted in Eig. 6. 
The diagrams show a very common situation where a primary 
electron excited by the laser pulse is loosing its energy on the 
way to the detector by exciting a secondary electron. There 
could be either bare or screened Coulomb interaction between 
the two electrons. In the latter case some resonant phenom¬ 
ena related to the excitation of e. g. plasmon are expected. 
The SPE case (Eigs. 6 (a,b)) is identical to the process of sec¬ 
ondary electron excitation considered by Caroli et al. [21]. 
All DPE processes covered by the diagram at Eig. 6 (c) form 
a subset of the SPE process. The only difference between the 
two scenarios is whether primary, secondary or both electrons 
are observed in the detector. It is obvious that one reduces 
the DPE diagram to the SPE ones by integration over the en¬ 
ergy and momentum of the secondary, or primary electrons, 
respectively. 

Since we do not take into account the interaction between 
the two emitted electrons (as given, for, instance by two F- 
blocks at Fig. 5 (a) one can express the final result for the cur¬ 
rent as a matrix element over the direct product of two single¬ 
particle scattering states. This is typically a good approxima¬ 
tion for the case when two electrons have different energies 
(momenta), or for approximately equal kj and k 2 in the case 
of larger energies [48]. 

To work this out consider a part of the DPE diagram 
that contains a product of two GEs involving the external 
momentum k. Introducing the Fourier representations for 
""5 = n G>,(-t') = 

f-cc expressing the interacting GF as a product 

of the Mpller operator and the free-particle Green’s function 
(see Appendix B) we obtain expressions similar to Eqs. (60). 
Thus, in the time domain the product of two interacting single¬ 
particle GFs reduces to a simple propagator computed on the 
scattering states with incoming boundary conditions: 

X 6»(-T)6i(-T')e^'’''"^'^ (68) 

As an exercise let us evaluate the diagram at Fig. 7 (a) de¬ 
scribing the SPE process with extrinsic plasmon losses. The 
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current is given by the following expression in the time do¬ 
main: 


/k = lim277lim 

rj—^Q 6~^0 






mUr, t') 

xG^,,(0, T')G>k(T, 0)G:;(f, T)(Aah)\ (69) 


Representing the lesser Green’s function on the vertical track 
in terms of the electron spectral function (normalized as 
Yjb JToo ^^bbiO - N, N is the number of electrons in the 
system) 


0 = ‘ £ ^A,b(Oe-‘^^''~'\ (70) 

and the greater component of the screened interaction in terms 
of the plasmon spectral function 

W>,,(T, r') = -i £ (71) 


representing time-ordered G^r^(T',t') and anti-time-ordered 
G^^(f, r) as Fourier integrals and using expression (68) we 
obtain: 


f‘ 


abed 

X \ d{cj\a)2)2ri 
1 


I 

1 


X 


(jj + ^ - oj\ - irj (jj + ^ - 0)2 + IT] 

-GAc((D2)GA(coi) 


oji - ^ - Sk - id OJ2 - ^ - Sk + id 

x{x^^V)AcdAdb(0(Aaby'ix\xV). 


(72) 


Now the limits can be taken making use of an identity discov¬ 
ered by C. O. Almbladh [16] (see Appendix E). It transforms 
the product of four fractions in the equation above into the 
product of three d-functions (2n)^d(oji - o) - 0(5(w2 - a> - 
+ fk - tu -1- ^), and after the frequency integration we 
obtain 


7k = '2.TT 


J-oc^TT Jo In 


d(^ + Sk-oj- A 


X J:i(xx')(Yl '^\x')BxA^){x\Xk ’> 
X [G“(w -h ^)AA(C)A^G++(w -h 


(73) 


The two-particle current is obtained along the same lines us¬ 
ing the energy flow as shown on Fig. 7 (b). 




"73 1 

8 

J-oc 27r ^ 

1-00 271 , 

o 

1 


¥ Ski -0)-A 

£{xx'zz')(xty)WAAm::i^){x\xt^) 


X (xL^\z}Az'z(0(z\xi^^}d(sk,_ 


X [G“(w -h OAA{OA^G^^{co + (74) 



FIG. 7. Energy flows in (a) SPE diagram with external plasmonic 
losses, (b) DPE diagram describing a related plasmon assisted pro¬ 
cess. Analytical expressions corresponding to these diagrams are 
first written in the time domain, then the integrations are performed 
by Fourier transforming all the propagators, and lastly the limits 
7/ —> 0 and d —> 0 are taken. 


Similarly to the previous case, the limits t; —> 0, b —» 0 yield a 
product (of five) 5-function which were subsequently used to 
perform three frequency integrations here (see Appendix E). 
All the quantities in Eqs. (73,74) can be expressed in terms 
of the spectral functions. We can, for instance, start with a 
general expresion for the time-ordered function in terms of 
functions on the Keldysh contour: 

r-(r) = fd(r) + 6(T)t(T) + 0(-r)fAT), (75a) 

t^(T) = -/5(t) + e(-T)f>(T) + 0(t)/<(t), (75b) 

where in the first equation t = t- - t'_ is equal to the time- 
difference on the forward branch of the contour, and t = 
is equal to the time-difference on the backward branch of the 
contour in the second equation. After the Fourier transform 
fioS) - f^dTe“^'^f(T), we have 


/--(m) = f + 


oo 


rdoj' 

irw) 

irw) 

j 2n 

OJ - oj' + id 

OJ — oj' — id 


—oo 


(76) 


The fluctuation-dissipation theorem at zero temperature al¬ 
lows to express the lesser and greater propagators in 
terms of the corresponding spectral functions (Kubo-Martin- 
Schwinger (KMS) conditions [10]): 

G^(«) = i0(fd - u)A{o)), G^(w) = -i0(to - ]i)A{o)). 


The screened interaction obeys KMS conditions for bosonic 
propagators: 

W^(«) = i0(-u)B{o}), W^(m) = -i0{oj)B(oj), 

with the symmetry property for the spectral function Bi—oS) - 
-B(to) (follows e. g. from the fact that W^(t, f) is a real func¬ 
tion or, more precisely, a Hermitian matrix). We have already 
used these equations (cf Eqs. (70,71)) to express SPE current 
in terms of spectral functions. Using Eq. (76) we can write the 
spectral representation of the fermionic propagator 


G-(oj) 



-7^A{o) ) 
2n 


0(ji - oj') 
OJ- oj' - id 


¥ 


0(oj' - jj) 
OJ — Oj' + id 
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where ji is the Fermi energy. The anti-time-ordered GF is ob¬ 
tained similarly = -[G (w)] ^ The screened interac¬ 

tion is expressed as an integral over the positive frequencies: 


W (io) = V -H 



dai' 

271 


B{jy 


0)^ 


2a)' 

- (co' - isy ’ 


while W^^iaS) - 

Let us consider plasmon-mediated DPE. This process is of 
relevance for metallic and large molecular systems. Since 
plasmon is a long wavelength or small momentum electronic 
excitation it is useful to go from the abstract basis to momen¬ 
tum representation and write W {k, co) in a short form as 


W(k, oj) - Vk 


1 H- 



(77) 


where u^(k) is the plasmon dispersion, cOp = Wp(0) is the clas¬ 
sical plasmon frequency, and is the matrix element of 

Coulomb interaction. It is clear that in this form the plasmon 
peak completely exhausts the /-sum rule. Such plasmon pole 
approximation for the screened interaction is broadly used in 
the electronic structure calculation when full-fledged calcula¬ 
tions are not feasible. Similarly, it can be used to simplify 
Eq. (74). 


C. Numerical results 


Let us make some simplifications. Usually it is a good ap¬ 
proximation to start with the mean-field Green’s functions 


g;;(w) - 

aeocc 


{x\a}na{a\y} 

to - Sa — id 


z 


{x\a}na{a\y} 

CO - Sa + id ’ 


(78) 


where ria is the occupation number of the state a and fia = 
1 - Ha- After straightforward, but tedious calculation the fre¬ 
quency integrations in Eq. (74) can be performed (for techni¬ 
cal reasons it is better to start from the time rather then fre¬ 
quency expression, and it can be obtained by directly tran¬ 
scribing the diagram at Eig. 7 (b) using standard rules) yield¬ 
ing the following expression for the two-particle current: 


Aikj = 47r V 


"I" ^b ^ki ^k2) 


^ki ^k2)(^ki ^k2 


z 

qiq4 


C(/k2rf] 

Cvq.ol 

f fid 

te) 

1 

Xsd - ekj)^ 

- Wp(?l) 

[(sd - ek2)^ 

- Wp(?2). 


• Sd) 

(79) 


with the following matrix elements 




(80) 


Notice that it is not necessary to separately treat the bare 
Coulomb interaction, it can be recovered as Wp —> oo limit 
as explained in [62]. 

Let us compare Eq. (79) with the general result obtained 
using the Eeshbach projection formalism (64). Eor the mean- 
field approximation (78) the two-particle spectral function is 


diagonal and is given by the convolution of two single-particle 
spectral densities: 

Jd'^Atb(^-OAdd(0 

= rihriddi^ - / - Eb)d(y - Ed) 

^nhtiddiEb+Ed-0- ( 81 ) 

The energy conservation for the whole process, which is given 
by the (5-function in the numerator of (79), is expressed in 
terms of the two-particle spectral function + ®k 2 “ w), 

(cf. Eq. (65)). The denominator of the first line reflects the 
resonant character of the considered two-step process. Erom 
the resonance conditions (zeroes of the denominator) we see 
that the double photoemission is enhanced when a and c are 
continuum states and therefore we denote them as and k^. 
We replace the scattering states and entering the 
matrix elements (80) by the plane-waves and perform the in¬ 
tegration yielding = (5(k-ka-q). Combining all together 
we obtain the following concise expression for the plasmon- 
assisted DPE process: 

7kik2 = ^TT A^^j(ek, + ek2 - w) 

k„k, bd 

<ki -H k2 - K\d}{d\iii -H k2 - kc) 

X - 

(Sk^ + Sd — £ki ~ fik2)(£ki + £k2 ~ Sk„ — Ed) 

X W(ki - kc, Ed - ek2)l^(ki - k^, Ed - ek2). (82) 

We have seen that the plane-wave approximation for the scat¬ 
tering states (i. e. the Mpller operator is given by the identity 
operator) results in a great simplification for the two-particle 
current: it is given by a sum over two bound states (they cor¬ 
respond to two lesser propagators in the diagrammatic repre¬ 
sentation of this process) and by the two momentum integrals 
corresponding to the propagators of the secondary electron. In 
contrast, in the full-fledged calculations based on Eq. (79) the 
momenta of the secondary electron and the emitted electrons 
are not rigidly related. Therefore, in general two additional 
momentum integrations are required. This will be the subject 
of a forthcoming publication where this formalism is applied 
to a large molecular system. 

The DPE process described by Eq. (82) is suited to probe 
the plasmon dispersion and damping. Eirst, let us look at 
the classical plasmon that carries vanishing momentum and 
otherwise is strongly damped. This leads us to consider the 
case kfl ss kc ~ ki, and Ed - Ey,^ - cOp is the condition for 
the plasmon resonance. In this case the second line reduces 
to Kk2|(7)P/a)p, and is clearly off-resonance. The situation 
greatly changes if we allow for the plasmon to cany finite mo¬ 
mentum qc and consider a large momentum of the secondary 
electron k^ ss k^ ss ki > -/Eijj. Eor simplicity take a sym¬ 
metric situation when both screened interaction lines carry 
approximately the same energy and momentum and denote 
K ss 2 (ka -I- ki) ss 5 (ke -H ki) and q ss ka - ki ss kc - ki. In 
this case one achieves the resonant enhancement when 

ek„ - Ski ek< - Ski = 2(q ■ K) = Wp. 
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Thus for colinear k^, kj. and ki the probability for the 
plasmon-assisted emission of the secondary electron is en¬ 
hanced when K reaches the value of 

In order to illustrate the features arising due to the plasmon- 
assisted process in an experiment, we computed the current 
for a simple model system. To be concrete, we consider the 
basic jellium model for the C(,o molecule (treated as spheri¬ 
cally symmetric) [63, 64], which is known for its pronounced 
(dipolar) plasmon resonance at a>p ~ 22 eV. Inserting a 
smoothed box-like potential as approximation to the Kohn- 
Sham potential, we solved the Schrodinger equation for the 
120 orbitals required (240 electrons in total). This procedure 
yields the single-particle energies sj associated to the orbitals 
ipd{v), from which we can compute all quantities in Eq. (82). 
Because of the spherical symmetry, we can separate the ra¬ 
dial and the angular dependence, that is ^6^(r) = 

{Y(m{r) are the spherical harmonics) and only solve the ra¬ 
dial Schrodinger equation. For the optical matrix elements, 
we choose the length gauge and assume a linear polarization 
along the z axis (A = z). Since we are not interested in the ab¬ 
solute scale a prefactor proportional to the field strength will 
not be included. The matrix elements Akz, attain the form 


Aki = 47r ^ 


Sbeik) 


tm 

J f>CC 

di 

0 


dr r^Uh(r)je(kr), 


where jt denotes the spherical Bessel function. The co¬ 
efficients ^6 obtained from the standard Clebsch- 

Gordan algebra [65, 66]. Similarly, the Fourier-transformed 
orbitals (k|(f) = can be expressed in terms of the 

Bessel transformation; 0d(k) = An'ud{k)Yt^miik) with'Mrf(k) = 

poo 

Jy drrud(r)je^(kr). 

Next we transform the summation over ka and k^ into in¬ 
tegrations and substitute them by the integration over the mo¬ 
mentum transfer vectors qa,c = ki - k„_c- At this stage, no fur¬ 
ther simplification can be made, such that the six-dimensional 
integral has to be evaluated. However, it is reasonable to 
consider c as small, since the plasmon branch enters the 
particle-hole continuum for growing momentum, where it is 
strongly damped. Hence, we introduce the momentum cutoff 
^max and assume ki,k 2 » qmax- Thus, we approximate Ak^ = 

Ak,-q„fc ~ Ak,fe and 0^(ki-i-k 2 -k„_£.) = (l>di ^2 + f\a,c) ~ ^^(kz)- 
Furthermore, we integrate over the spherical angles of ki and 
k 2 , keeping only the dependence on their magnitude. Thus, 
the two-electron current can be written as 


EE [^^(^ 2 )! 


bd km 


1 -I- Re 




(Sd - Sk2 - iY)^ 


(83) 


Fd{ki,k2), 


where 



dq 


1 

q^ + 2k\q - k^+ 2ed 





FIG. 8. (Color online) The symmetrized two-electron current as a 
function of the photoelectron energies (energy-sharing diagram) for 
typical parameters: tu = 2.0 and tjp = 0.8. The color scale is the 
same for all three panels and runs from dark blue to red, indicating in¬ 
creasing values, (a) The process is mediated hy the pure Coulomb in¬ 
teraction. (b) Pure plasmonic contribution, (c) Total (bare Coulomb 
and plasmonic contributions) signal including the interference terms, 
(d) Equal energy sharing (si;, = ek 2 ) for the current and trace of the 
two-particle spectral density (shaded curve). 


Note that we inserted the imaginary shift /F in the energy ar¬ 
gument accounting for a finite width (lifetime in the time do¬ 
main) of the plasmon resonance (which is assumed disper¬ 
sionless for simplicity). 

In an experiment, the distinction between primary (ki) and 
secondary electron (k 2 ) is, of course, not possible. For this 
reason, the photo-current needs to be symmetrized (let us de¬ 
note it by Representing the 7*^™ as a function of Sk, 

and £k 2 yields the typical energy-sharing diagrams (Fig. 8). 
Spectral properties of the system (dominated by AY-^s)) dis¬ 
play themselves along the main diagonal, as only the sum 
Ski + ek 2 enters. Dominant scattering events mediated by the 
(screened) interaction on the other hand are visible along lines 
fiki = const, (or ek 2 = const.). As Eq. (83), indicates the 
two-particle current contains contributions from (i) the bare 
Coulomb (two interacting lines in Fig. 7 (b) are not screened), 
(ii) plasmonic scattering (both lines are screened), and the in¬ 
terference terms, (i), (ii) and the total contribution is shown at 
panels (a), (b) and (c) of Fig. 8, respectively. For vanishing T 
the current is dominated by sharp plasmonic resonances. For 
finite damping parameter such as used for the present simula¬ 
tions (T = 0.1, we use a realistic value as in Ref. [67]) the in¬ 
terference terms are important; we still have a large plasmonic 
contribution (viz. Fig. 7 (b)), however, the bare Coulomb con¬ 
tributes with the opposite sign. Therefore, in total current the 
large peak at ek 2 ~ 0.15 becomes less pronounced and ad- 


Fdikukt) = 
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ditional peaks at higher energies (e.g. at ~ 0.5) appear. 
The whole spectral width of the signal is limited by the two- 
particle spectral function shown at Fig. 7 (d) as a shade curve. 
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VI. CONCLUSIONS 


There is a large number of theoretical works devoted to 
the interaction of light and matter which involves the emis¬ 
sion of one or more electrons. This contribution is meant to 
expose parallels between the single and the double electron 
photoemission in a formal way. We started by defining cor¬ 
responding observables and deriving expressions for one- and 
two- particle currents based on the first-order time-dependent 
perturbation theory. These expressions are suitable if exact 
formulas in terms of many-body states are required. In order 
to obtain computationally useful expressions many-body ef¬ 
fects should also be accounted for in a perturbative fashion. 
Thus, in the first part of the manuscript we applied the pro¬ 
jection operator formalism. Starting from the explicit form of 
the projection operators dividing the whole Hilbert space of 
the system into that of the emitted electron(s) and the target 
we derived the effective one- and two-particle Hamiltonian, 
discussed integral equations for the Green’s functions describ¬ 
ing emitted particles and demonstrated a close connection of 
this formalism to the nonequilibrium Green’s function theory. 
For the latter, one can easily derive the diagrammatic expan¬ 
sions for one- and two-particle currents starting from the time- 
dependent perturbation theory and using the adiabatic switch¬ 
ing of the electron-electron interaction. Hence, we have elec¬ 
tromagnetic field switched on at the remote past (as and 
independently adiabatically switched on the interaction such 
that the total Hamiltonian takes a form Hg - Hq + 

We analyzed in details the diagrammatic structure of one- and 
two-particle currents. It is surprisingly simple: one starts 
with the density-density response function which neces¬ 
sarily contains two blocks associated with the forward 
and backward (“H-”) parts of the Keldysh contour. Request¬ 
ing that one or two lines flowing from to “+” blocks are 
associated with scattering states (with momenta k,) one ob¬ 
tains exactly the diagrams for SPE and DPE currents showing 
the close connection between these types of light-matter in¬ 
teraction. It is not difficult to generalize this approach to an 
arbitrary number of particles. Einally, we presented a detailed 
analysis of the plasmon-assisted DPE and showed that if one 
of the emitted particles is unobserved, its diagrammatic repre¬ 
sentation reduces to the one describing external losses in the 
SPE process considered by Caroli et al. [21]. Plasmon pole 
approximation was employed to derive computationally man¬ 
ageable expressions. We illustrated the distinct features to be 
expected in an experiment by analyzing the simple and yet re¬ 
alistic jellium model for the Ceo molecule. This will extended 
used in the forthcoming paper devoted to the ab-initio treat¬ 
ment of this large molecular system. 


Appendix A: Particle-impact ionization 


Under some circumstances the formalism developed in the 
main text can be extended to other mechanisms of ionization, 
e.g., particle-impact ionization. The basic requirement we im¬ 
pose is the distinguishability of the projectile from the target 
electrons. This applies also for a projectile electron if the im¬ 
pact energy is high and the small momentum transfer is small 
(optical limit). 

The target we describe by the Hamiltonian Eq. (4). The 
Coulomb-interaction between the projectile (with charge Z) 
and the sample reads 

^ ~ 2 ^ i ^ 1 ’^avbfiCadld^Ch . (Al) 

ab vfi 

dy (dl) is the annihilation (creation) operator of the projectile 
states |v). These states can be chosen as the eigenstates of the 
projectile Hamiltonian with energy Sy. 

Assuming that the projectile initially possesses the momen¬ 
tum k;, we can construct the asymptotic state prior to the in¬ 
teraction (that is, at f = -oo) as the product state 

I'Eo.k,) = |k,) ® it'd) . 


V is switched on reaching its full strength at f = 0. As¬ 
suming that this strength is still way smaller than the inter¬ 
nal interaction within the sample, we can apply the first order 
perturbation theory (i.e., the first Bom approximation in the 
projectile-target interaction [47]). Denoting the full Hamilto¬ 
nian hy H hp one may write 

= I'I'o.k.) + lim--VIT'o.k,) . (A2) 

Eo -H fik, - H - hp-\- irj 

The projectile has a well defined final momentum k/. In anal¬ 
ogy to Sec. IV A we introduce the particle number operator 

Ak ^ P/4E/ 


with Pf - |ky^)(ky| projecting only onto the projectile space. 
Ak acts on the system’s states only (including the ejected elec¬ 
trons upon particle impact). Evaluating then the current as 
in Sec. IIA and approximating the projectile states by plane 
waves (r|k) = yields 


7k = lim277(T'o,k, 


1 


Eq + Sk, — H - hp — irj 


-P/CkCkP/ 


Eq -h fik, - El - hp + irj 

lim277(T'o|y‘=®(q) 


r^'I'o.k,) 


1 


(A3) 


1 


Eq -H fik, - fik, - 77 H- irj 


Eq + fik, - fik/ -H-irj 


-y^‘f(q) T'o) , 
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where q = k, - k/ is the momentum transfer, and y‘^®(q) is 
the effective single-particle operator acting on the target, ex¬ 
plicitly 

y^ff(k,- - kf) = <k,|y|ky) = I ^ VaK.hkdcb ■ (A4) 

ah 


In this optical limit 


y 




(A5) 


k. two-particle case: For DPE the two-particle Green’s 
function over the excited state 'Vjf is required; 


--(pp) 

^pq, 




1 


z-Hp- Sp(z) 


-pt J 




where the projection operator is defined by Eq. (35). This 
propagator can be related to the scattering Green’s function of 
the two-particle system in the presence of the optical potential 
of doubly-ionized target: 


acts similar to the light-matter interaction A; the transferred 
energy (or energy loss) et; - Sk^ resembles the photon energy. 


Appendix B: Green’s functions 


with e® - Eq - ^pq^kik 2,8 likewise expressed in 

the form (Bl). 


Let us recast the following many-body correlators from Appendix C: Matrix identities 

Sec. IV A 


The formalism of presented here works in finite as well as 
in infinite-dimensional Hilbert spaces. Eor illustration we for¬ 
mulate it in the matrix form. Given M is square block matrix: 

in the form of one-particle averages. We define the particle 
propagator of one-particle system in the presence of optical 
potential W(z)'. 


"B 

C £) 


(Cl) 




1 


• Hp - E„(z) 




Consider Gjfq „(w -i- Eq ± if)- The matrix element of the effec¬ 
tive Hamiltonian operator in its definition can be simplified 
to 

<'I':|cp(i/p + Ep(z))4|T':) = e: + {p\Hf + Wa(z)\q}, 

where we decompose the total V-particle Hamiltonian // as a 
sum of three terms: 

H^Hf + H^ + V. 

Here Hj is the free-particle Hamiltonian, is the Hamilto¬ 
nian of ionized system: 

and y is the frequency independent part of the self-energy. If 
the optical potential is identified with the self-energy then one 
we can relate two propagators 

Cpq,Q'(tU -t Eq it irj) — ^pq^Q,(tU -t Sa), 

where we introduced the Green’s functions ^pq\w) = 

S^pqio) ± irj) and Sa - Eq - E^. Erom the formal scatter¬ 
ing theory (see Sec. 20 of Joachain [47]) and independent of 
the concrete choice of the representation we can express them 
in terms of the Mpller operator and the free-particle Green’s 
function 

(Bl) 


where D is square invertible matrix the Schur comple¬ 
ment [ 68 ] (also known in physics as the Eeshbach map [57, 
69, 70]) is defined as: 

5^ = - SD ‘C. 

We might think of A1 as a Hamiltonian operator acting in 
some larger Hilbert space, whereas JTI is the same operator, 
but acting in a physically relevant subspace. Be P the projec¬ 
tion operator onto this subspace (PAiP - J?[) and Q - I — P 
is its complement {QMQ - D). Eor definiteness we may 
take Af to be a compact self-adjoint operator on the Hilbert 
space describing an V-fermion system and JTI its projec¬ 
tion upon the Hilbert space of two particles Because of 
the couplings between subspaces (for physical Hamiltonians 
obviously holds B = C' ) A1 and JTI have different spectral 
properties. Nonetheless, one can show the following equiva¬ 
lence 

Afy = 0 ^PV = 0 , (C 2 ) 

for a vector V e If A1 = H - El the first part implies 

that y is an eigenvector of H with the energy E. The second 
part implies that Py is a corresponding eigenvector of J^(E) 
with the same energy: 

{Hp H- Sp(E) - EIp)PV = 0. (C3) 

Expression for the self-energy (27) is derived for instance 
in Sec. 20.2.3 of Joachain [47]. A mathematically rigorous 
proof of the theorem (C2) as well as other properties of the 
Eeshbach-Schur map can be found in Chap. 11 of Gustafson 


g\co) = 
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and Sigal [71]. It is further possible to write the inverse of the 
matrix Af explicitly [72]: 

M-i - f 

This identity is natural to apply to compute resolvents. For in¬ 
stance, Eq. (26) is given the first line of Eq. (C4). This formula 
can also be found in Almbladh as Eq. (19) [16]. 


Using this result and the vacuum assumption for the initial 
states we can compute a matrix element entering the Fermi 
golden rule formula for SPE: 

{'P^IckAl'Po) = 2 A,i<'P:|ckclci|'I'o> 

ab 

^ ^(klAlfoX'T^lcil'Po) - <k|A|,^,>, (D7) 

b 

and DPE: 


Appendix D: Properties of projection operators 

Our basic assumptions for operators with continuum in¬ 
dices CpI'P^) = 0 and Cp|'P^^) = 0 imply that final states of 
the target are the vacuum states for these operators. Thus, 
standard Wick’s theorem can be used for the calculation of 
various correlators. It follows 

(Dl) 

Ck,Ck,44|'PjX = (4,p4.q - 4.q5k.p)|T|^)- (D2) 

These equations lead to the indempontency relations PaPa - 
Pa and P/sPfi = PfS and to the properties 

cIl'i'aX'i'.Ick = P.4ckP„, (D3) 

= Ppclclw^Pfi■ (D4) 

The matrix element of a one-particle operator O = 0{x\) + 
0 {x 2 ) over the determinant two-particle states {xxX 2 \ab) — 
-^{(pa{x\)4’b{x2)-(pb{xi)4>a{x2)) Can be verified by direct eval¬ 
uation: 


{ab\d\cd) = {a\d\c)5bd + {b\d\d)6ac 

- {a\d\d)6hc - {b\d\c)5ad. (D5) 

If one of the states is a two-hole Dyson orbital the matrix ele¬ 
ment is computed similarly: 


<'p2+|ck,Ck,A|T'o> - ^ A^X'P^Ick.Ck.clcIT'o) 

ab 

^ ^ [<ki|A|f7)4,. - <k2|A|f7)(5k,.]<'i'^1ciQ|'Po> 

be 

= <kik2|A|4">>. (D8) 

We used an assumption Ckl'Po) ss 0 to derive (D7) and 
CkiCkjl'Po) ~ 0 to derive (D8). 

Appendix E: Sokhotski-Plemelj-type identities 


Eollowing identities were used to perform frequency inte¬ 
grations leading to Eqs. (73,74). 


1 


1 


1 


lim lim Irj- 

tj-^os^Q a>i - zi - irj ci>2 — Z2 + irj — id 

3 


X- 


1 


Z3 + Z2 - Zi - + id 

for co\ - a>2, and 


Y]27rd(zi-cod, (El) 


1=1 


1 


1 


1 


lim lim Itj- 

Tj-tO 6^0 coi - z\ - iri 0)2- Z2 + irj Z3 - Z2 + bL>2- 0J3 - id 

1 1 

X- 

Z4 - Zl + 0)1 - 0)4 + id 0)4 + 0)5 - Z4 - Z5 + id 

1 


- - = n 2nd(zi - oJi), (E2) 

— 7 n — — in 11 


0)3 + 0)5 - Z3 - Z5 - id 


{ab\d\^f^) = i J]{ab\d\cd}{^'l^\CaCd\^o} 

cd 

= ^ {{a\6\c)dbd - {b\6\c)dad){%^\ccCd\'^^{). (D6) 

cd 


for 0)1 - 0)2, 0)3 - 0)4. The first equation appears in [16]. To 
the best of our knowledge the second equation has not been 
addressed in the literature. These identities can be verified by 
the Eourier transformation with respect to z, variables. 
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